home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 107_01 / clogs.c < prev    next >
Text File  |  1984-06-03  |  18KB  |  652 lines
  1.  
  2.  
  3.   /*  CLOGS.C       *****
  4. A group of programs in C, using the BDS-C floating point package
  5. as modified by LCC (FLOAT+44.*) and depending on the ability to
  6. insert null characters in a string available in BDS-C V 1.44.
  7. Four functions are handled:
  8.     log10, exp10, expe, pi
  9. In addition, there is a service function exprange() which returns
  10. a false (00) if the exponent of a floating point variable is reaching
  11. near the end of the usable range.
  12. These are discussed in detail in CLOGS.DOC
  13.  
  14. L. C. Calhoun
  15. 257 South Broadway
  16. Lebanon, Ohio 45036  513 day 433 7510 nite 932 4541
  17.  
  18.   */
  19.  
  20.  
  21. char *pi(result)
  22.  
  23.  /* result is a standard character array char result[5]; in calling
  24. program as used for floating point variables.  The return is a pointer
  25. to result, and the value of pi stored in the result array in floating
  26. point */
  27.  
  28. char *result;
  29. {
  30.  char *piconst, *fpasg();
  31.  
  32.  piconst = "\171\356\207\144\2";
  33.  return (fpasg(result,piconst) );
  34. }
  35.  
  36. char *expe(result,xin)
  37.  
  38.  /* computes the base of the natural log "e" raised to the "x'th"
  39.     power.  Checks are made for out of range values and result is
  40.     defaulted to 0, 1, or a large number as appropriate.  There are
  41.     no error flags.  The arguments are floating point character
  42.     arrays char result[5], x[5]; in calling program.  Return is
  43.     a pointer to result, and "e to the x" stored in result.
  44.  */
  45.  
  46. char *result, *xin;
  47. {
  48.     char *zero, *one, *large, *coef[7], *eghty6, *meghty6;
  49.     char intres[5], xint[5], x[5];
  50.     char  *fpmult(), *fpasg(), *fpdiv(), *fpchs();
  51.     int signx, index, bigexp, smallexp, zeroin;
  52.     int exprange();
  53.  
  54.     bigexp = smallexp = zeroin = 0;
  55.  
  56.     zero ="\0\0\0\0\0";
  57.     one = "\0\0\0\100\1";
  58.     large = "\0\0\0\100\175";
  59.     eghty6 = "\0\0\0\126\7";
  60.     meghty6 = "\0\0\0\252\7";
  61.  
  62.     coef[0] = "\0\0\0\100\1";
  63.     coef[1] = "\140\326\377\177\376";
  64.     coef[2] = "\130\373\3\100\374";
  65.     coef[3] = "\200\1\352\124\370";
  66.     coef[4] = "\351\253\362\131\364";
  67.     coef[5] = "\21\213\32\133\357";
  68.     coef[6] = "\371\330\260\134\354";
  69.  
  70.   /* preserve input datum */
  71.     fpasg(x,xin);
  72.  
  73.  /* check for sign */
  74.     if (x[3] < 128)   /* positive */
  75.        {signx = 1;
  76.          }
  77.     else
  78.          {signx = 0;
  79.           fpchs(x,x);
  80.          }
  81.  
  82.    /* check for zero and out of range of fp var */
  83.  
  84.    /* check for zero and very small numbers */
  85.     if ( ((x[4]>127) && (x[4]<226)) || ( (x[4]==0) &&
  86.         (x[3]==0) && (x[2]==0) && (x[1]==0) && (x[0]==0) ) )
  87.         zeroin = 1;
  88.  
  89.    /* check for very small exponent */
  90.     if ( fpcomp(xin,meghty6) < 0)
  91.         smallexp = 1;
  92.  
  93.   /*  check for very large exponent */
  94.     if ( fpcomp(x,eghty6) > 0 )
  95.         bigexp = 1;
  96.  
  97.   /*  now if small number or zero, result is one */
  98.   /*  now if big number and positive, result is large number */
  99.   /*  now if big number and negative, result is zero */
  100.  
  101.     if (zeroin) return (fpasg(result,one) );
  102.     if (smallexp) return (fpasg(result,zero) );
  103.     if (bigexp) return (fpasg(result,large) );
  104.  
  105.  /*  all exceptions taken care of, so evaluate rest  */
  106.         fpasg(result,one);
  107.         fpasg(xint,x);
  108.         index = 1;
  109.         while ( (index<=6) && exprange(xint) )
  110.            {
  111.            fpmult(intres,coef[index],xint);
  112.         fpadd(result,result,intres);
  113.         fpmult(xint,xint,x);
  114.         index++;
  115.            }
  116.   /* now do the square square */
  117.     fpmult(result,result,result);
  118.     fpmult(result,result,result);
  119.  
  120.   /* now treat sign appropriately */
  121.     if (signx) return (result);
  122.     else
  123.        {fpdiv(result,one,result);
  124.         return (result);
  125.        }
  126. }
  127.  
  128. char *exp10(result,xin)
  129.  
  130. /* similar to expe, except result returned is 10 raised to the x
  131. power    the antilogarithm to base 10  */
  132.  
  133. char *result, *xin;
  134. {
  135.     char *zero, *ten, *one, *large, *thty8;
  136.     char xint[5], *coef[7], intres[5], tenfac[5], x[5];
  137.     int index, bigexp, smallexp, signx, tenpower;
  138.     int exprange();
  139.  
  140.     bigexp = smallexp = 0;
  141.  
  142.     zero ="\0\0\0\0\0";
  143.     one = "\0\0\0\100\1";
  144.     large = "\0\0\0\100\175";
  145.     ten = "\0\0\0\120\4";
  146.     thty8 = "\0\0\0\114\6";
  147.  
  148.     coef[1] = "\65\264\256\111\1";
  149.     coef[2] = "\0\14\330\124\0";
  150.     coef[3] = "\0\46\354\100\377";
  151.     coef[4] = "\24\140\107\115\375";
  152.     coef[5] = "\242\304\361\155\372";
  153.     coef[6] = "\361\143\246\134\371";
  154.  
  155.  /* preserve input datum */
  156.     fpasg(x,xin);
  157.  
  158.  /* check for sign */
  159.     if (x[3] < 128)   /* positive */
  160.        signx = 1;
  161.     else     /* negative */
  162.        {signx = 0;
  163.         fpchs(x,x);
  164.        }
  165.  
  166.   /* check for very small or large numbers, check by exponent size */
  167.   /* check for zero or small */
  168.     if ( ((x[4]>127) && (x[4]<226)) || ( (x[4]==0) &&
  169.         (x[3]==0) && (x[2]==0) && (x[1]==0) && (x[0]==0) ) )
  170.         smallexp = 1;
  171.    /* check for big number */
  172.     if ( fpcomp(x,thty8) > 0)
  173.         bigexp = 1;
  174.  
  175.  /* if value is small or zero, return 1 as with expe */
  176.  /* if value is large and positive, return a large number */
  177.  /* if value is large and negative, return zero */
  178.  
  179.     if (smallexp) return (fpasg(result,one) );
  180.  
  181.     if(bigexp && signx) return (fpasg(result,large) );
  182.  
  183.     if(bigexp && !signx) return (fpasg(result,zero) );
  184.  
  185.  /* now reduce range of x to between zero and one */
  186.  
  187.     tenpower = ftoit(x);
  188.     itof(tenfac,tenpower);
  189.     fpsub(x,x,tenfac);
  190.     fpasg(tenfac,one);
  191.     while (tenpower)
  192.         {fpmult(tenfac,tenfac,ten);
  193.          tenpower--;
  194.         }
  195.  /* now evaluate series  */
  196.     fpasg(result,one);
  197.     fpasg(xint,x);
  198.     index = 1;
  199.  
  200.         while ( (index <= 6) && exprange(xint) )
  201.            {fpmult(intres,coef[index],xint);
  202.         fpadd(result,result,intres);
  203.         fpmult(xint,xint,x);
  204.         index += 1;
  205.            }
  206.  
  207.  /* now square result (note error in referenced article) */
  208.     fpmult(result,result,result);
  209.  
  210.  /* now check sign and make proper return */
  211.     fpmult(result,result,tenfac);  /* scale back up */
  212.  
  213.     if (signx) return (result);
  214.  
  215.     else return ( fpdiv(result,one,result) );
  216. }
  217.  
  218.  
  219. char *log10(result,sign,xin)
  220.  
  221.  /* computes briggsian logarithm of x which is a char[5]
  222.     floating point number.  Return is logarithm in result[5],
  223.     and sign pointed to by sign.  The logarithm is taken
  224.     of the magnitude, and sign information preserved
  225.     as required by sign.
  226.  */
  227. char *result, *xin;
  228. int *sign;
  229. {
  230.     char *zero, *ten, *one, *large;
  231.     char *sqrtten, x[5];
  232.     char gamma[5], gamnum[5], gamden[5], *coef[5];
  233.     char *half;
  234.     char intres[5], gamint[5];
  235.     int tenpower;
  236.     int index, bigexp, smallexp, signx;
  237.     int exprange();
  238.  
  239.     bigexp = smallexp = 0;
  240.  
  241.     zero ="\0\0\0\0\0";
  242.     one = "\0\0\0\100\1";
  243.     large = "\0\0\0\114\6";
  244.     ten = "\0\0\0\120\4";
  245.     half = "\0\0\0\100\0";
  246.     sqrtten = "\304\145\61\145\2";
  247.  
  248.     coef[0] = "\362\6\56\157\0";
  249.     coef[1] = "\30\346\21\112\377";
  250.     coef[2] = "\100\55\344\132\376";
  251.     coef[3] = "\106\73\244\140\375";
  252.     coef[4] = "\174\5\367\141\376";
  253.  
  254.  /*  preserve input variable */
  255.     fpasg(x,xin);
  256.  
  257.  /* check for sign */
  258.     if (x[3] < 128)   /* positive */
  259.        signx = 1;
  260.     else     /* negative */
  261.        {signx = -1;
  262.         fpchs(x,x);
  263.        }
  264.  
  265.   /* check for very small or large numbers, check by exponent size */
  266.   /* check for zero or small */
  267.     if ( ((x[4]>127) && (x[4]<209)) || ( (x[4]==0) &&
  268.         (x[3]==0) && (x[2]==0) && (x[1]==0) && (x[0]==0) ) )
  269.         smallexp = 1;
  270.    /* check for big number */
  271.     if ( (x[4]  >47) && (x[4] < 128) )
  272.         bigexp = 1;
  273.  
  274.  /* if very small, return - a large number
  275.     if very large, return + a large number  */
  276.  
  277.     if (smallexp)
  278.         {*sign = signx;
  279.          return (fpchs(result,large) );
  280.         }
  281.     if (bigexp)
  282.         {*sign = signx;
  283.          return (fpasg(result,large) );
  284.         }
  285.   /*  now bring into range 1 <= x < 10 */
  286.     tenpower = 0;
  287.     while ( fpcomp(x,ten) >= 0)
  288.            {tenpower++;
  289.         fpdiv(x,x,ten);
  290.            }
  291.     while ( fpcomp(x,one) < 0)
  292.            {tenpower--;
  293.         fpmult(x,x,ten);
  294.            }
  295.  
  296.   /* now if exactly one, no need to evaluate */
  297.     fpsub(gamnum,x,one);
  298.     if (((gamnum[4]>127)&&(gamnum[4]<209)) || ((gamnum[0]==0) &&
  299.        (gamnum[1]==0) && (gamnum[2] == 0) && (gamnum[3] == 0) ) )
  300.            {*sign = signx;
  301.         itof(result,tenpower);
  302.         return (result);
  303.            }
  304.  /* now compute gamma  for series  */
  305.  
  306.     fpsub(gamnum,x,sqrtten);
  307.   /* now check for size of numerator */
  308.     if (((gamnum[4]>127)&&(gamnum[4]<209)) || ((gamnum[0]==0) &&
  309.        (gamnum[1]==0) && (gamnum[2] == 0) && (gamnum[3] == 0) ) )
  310.            {itof(result,